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Summary. The moment method is a well known mode identification technique in asteroseis- 
mology (where 'mode' is to be understood in an astronomical rather than in a statistical sense), 
which uses a time series of the first 3 moments of a spectral line to estimate the discrete os- 
cillation mode parameters e and m. The method, contrary to many other mode identification 
techniques, also provides estimates of other important continuous parameters such as the in- 
clination angle a, and the rotational velocity Ve. We developed a statistical formalism for the 
moment method based on so-called generalized estimating equations (GEE). This formalism 
allows the estimation of the uncertainty of the continuous parameters taking into account that 
the different moments of a line profile are correlated and that the uncertainty of the observed 
moments also depends on the model parameters. Furthermore, we set up a procedure to take 
into account the mode uncertainty, i.e., the fact that often several modes {£, m) can adequately 
describe the data. We also introduce a new lack of fit function which works at least as well 
as a previous discriminant function, and which in addition allows us to identify the sign of the 
azimuthal order m. We applied our method to the star HD181558 using several numerical 
methods, from which we learned that numerically solving the estimating equations is an in- 
tensive task. We report on the numerical results, from which we gain insight in the statistical 
uncertainties of the physical parameters involved in the moment method. 

Keywords: Generalized estimating equations, time series, sandwich estimator, astrostatistics, 
discriminant function 

1. Introduction 

Stars consist of a number of gas layers with different temperatures, pressures, and chemical 
compositions. During their sojourn on the main-sequence, i.e., when they transform hydro- 
gen into helium, some stars are subject to oscillations which in turn provide astronomers 
with a wealth of information about the stellar interior. This is the subject of asteroseis- 
mology. Such oscillations typically exhibit multiple frequencies and manifest themselves at 
the surface of the star through variations in brightness, temperature, and surface velocity; 
some of these are observable. A star can oscillate in one or more of its "natural" frequencies 
determined by the internal structure of the star. With suitable inversion techniques it is 
possible to use the observed frequencies to derive information about this internal structure. 

To do so, however, the characteristics of the oscillations need to be considered first. That 
is, a mode identification (note that 'mode' is an astronomical term) has to be carried out, in 



2 J. De Ridder, G. Molenberghs and C. Aerts 



which one estimates the parameters characterising the oscillations from observational data. 
There are few mode identification techniques, and the properties of their estimators are 
rarely studied. Statistical uncertainties of the estimates, for example, are never reported. 
Nevertheless, from an astrophysical point of view, such uncertainties are important because 
wrong mode identifications can bias inversion techniques. It is therefore necessary to know 
a priori the extent of possible errors in the estimates. 

In this paper, we study the statistical properties of one particular mode identification 
technique, the so-called moment method. For examples of applications of this method we 
refer to, e.g., Aerts et al. (1998), Uytterhoeven et al. (2001), Aerts & Kaye (2001), and 
Chadid et al. (2001). 

2. Astrophysical Background 

As in any inferential method, the moment method uses a theoretical model to describe the 
observations. To understand the statistically relevant properties of this theoretical model 
and its parameters, we first briefly discuss some of the physics of stellar oscillations and 
how they are observed. 

Figure Ogives a diagrammatic illustration of an orthographic planar projection of the 
surface of an oscillating star, i.e., parallel with the line of sight. For our application, the 




Fig. 1 . A diagrammatic illustration of an orthographic planar projection of the surface of an oscillating 
stars for the mode {£, m) = (5, 3). In the left picture we look at the equator, and in the right picture 
we look almost at the pole of the star. 

most important aspect of stellar oscillation is the surface velocity. The lighter parts of the 
stellar surface have an inward velocity while the darker parts have an outward velocity. The 
figure is only a snapshot: the star varies periodically and half an oscillation cycle later the 
situation is reversed with the lighter parts moving outward and the darker parts moving 
inward. For slowly-rotating oscillating stars, each of the oscillation modes can be described 
with a spherical harmonic F/" (which is actually the basis of our illustration in Figure ^ , 
where £ is the total number of nodal lines and m is the number of nodal lines perpendicular 
to the equator. In reality, the motion of a surface element is more complex because it moves 
horizontally as well as vertically. 

In terms of model parameters we have to estimate 3 unknown parameters per oscillation 
mode: 2 discrete parameters and 1 continuous parameter. The 2 discrete parameters are 
the mode numbers £ and m of the spherical harmonic, which describe the configuration of 
the inward and outward going regions. To describe the 3-dimensional motion of the stellar 
matter, only one parameter is needed: the amplitude Vp of the vertical motion, since there is 
a theoretical linear relation between the amplitude of the vertical motion and the amplitude 
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Table 1. A summary of all unknown relevant model parameters, with their mean- 
ing and their physical range. The notation (3 = [Pi, 132, 133,13a)' for the continuous 
parameters will be introduced and used in Section|31 



Parameter 


Meaning 


Physical Range 


l 




Degree of the spherical harmonic 


{0,1,2,..-} 


m 




Azimuthal order of the spherical har- 
monic y/" 


{-^,---,o,---,+n 


Vp 


= l3i 


Velocity amplitude of the oscillation 


> 


a - 


= 132 


Width of line profile in absence of pul- 
sation and rotation (nuisance) 


> 


Ve 


= 03 


Equatorial rotational velocity 


> 


a : 


= I3a 


Inclination angle of the star 


[0°,360°) 



of the horizontal motion. To compute the constant of proportion K , however, the mass and 
the radius of the star are required and these quantities are often not very accurately known. 
Nevertheless, in what follows we will assume, as a first approximation, that K is known, to 
considerably simplify the treatment. 

A further continuous parameter related to the oscillation is the oscillation period P. 
However, for good datasets, this oscillation period can often be quite accurately determined 
from the data with other methods so that it is usually regarded as known. 

In the model, 2 additional unknown parameters not connected to the oscillations are 
present. A first one is the rotational velocity at the equator of the star, usually denoted 
as Ve- The second one is the inclination angle a under which we observe the star. This 
is illustrated in Figure ^ Both pictures show the same Y^™, but on the left hand side we 
are looking on the equator, while on the right hand side we are looking almost on the pole. 
Clearly, a has a large impact on how the surface velocity field is observed. 

A last unknown model parameter is specifically related to the kind of observational data 
we use. In the case of the moment method it concerns high-resolution spectroscopic data. 
The gathered star light is decomposed into its colours so that a detailed spectrum can be 
constructed, i.e., received light flux as a function of the wavelength of the light. At certain 
wavelengths, such a spectrum contains absorption lines where the light has been partially 
blocked by certain chemical elements at the surface of the star. An example of the Si"*" 
absorption line at A = 4f 2.805 nm for the non-radially oscillating star HD181558 is shown 
in the left hand panel of Figure|21 Here, an observational time series of TV = 30 high-quality 
spectra gathered by De Cat and Aerts (2002) is shown. The oscillations in the star cause 
the absorption line to change its position and shape in time. Precisely these line profile 
variations are used to estimate the parameters mentioned above. To model them, another 
unknown parameter is needed, denoted by a, which is related to the width the line profile 
would have in the absence of pulsation. From an astrophysical point of view this is an 
unimportant nuisance parameter. For convenience. Table ^summarizes all unknown model 
parameters mentioned above, their meaning and their physical range. 

Modeling the line profiles themselves turns out to be very computationally expensive. 
That is why Balona (1986) devised the moment method, which replaces each line profile by 
the first, second, and third moment denoted by j/i, 2/2, and 2/3 respectively. These quantities 
are measures for the average position, the square of the width and the skewness of the line 
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Fig. 2. In the left panel, a time series of Si+ (412.805 nm) absorption lines of the non-radially 
oscillating star HD181558 is shown. The line profiles are 'sorted' to cover an entire oscillation cycle 
of the dominant mode, which has a period of about 29h42m. Each of the line profiles is vertically 
shifted to obtain a clear visual effect. In the right panels, the first moment yi (in km/s), the second 
moment y2 (in km^/s^) and the third moment (in km^/s'') of all line profiles are shown as a function 
of the oscillation phase 0. 



profile. Precisely, 



/+^[l-p(</.,A)]A"dA 



(1) 



where p{(f>, A) is the line profile function at phase (j) and for wavelength A. For each (f>, there 
is a separate line profile in the left hand panel of Figure El leading to a point for each of the 
right hand side sub-panels, corresponding to n — 1,2,3. In practice, no higher moments 
are considered since these are often noisy and unduly complicate the calculations. One 
commonly expresses the moments in (km/s)", by transforming the wavelength A in |^ to 
a velocity using the Doppler transformation formula: 

A-Ao 

V = c , 

Ao 

where c is the speed of light. The moments yni'f) can be expressed in terms of time t 
as well, where is defined as imodP/P, with 'mod' standing for the decimal part and P 
the oscillation period. A time series of theoretical moments can be computed much faster 
than one of theoretical line profiles. The nuisance parameter cr, however, remains. In the 
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right hand panels of Figure |21 we show a time series of the three moments for the star 
HD181558. The computation of such moments from spectral line profiles take the form of 
intensity-weighted sums, sums of squares, and sums of cubes. 

Fifteen years after the introduction of the moment method, this mode identification 
technique is still very relevant (see the recent references given in Section Indeed, the 
effort required for direct computation of time series of line profiles is currently still too 
computationally demanding to be useful for mode identification, even under simplifying 
assumptions concerning the shape of the absorption line. 

A theoretical moment at one point of time is computed by integrating over the contri- 
butions of all points on the visible stellar surface. Closed-form expressions for the moments 
exist (Aerts et al., 1992), but they are quite lengthy and of little practical use to computing 
derivatives, such as in Sectional We opt for different, computationally more advantageous, 
expressions, which involve an integration with bounds depending on the inclination angle 
a. 



3. Current Statistical Status 

The moment method is a multi-response problem where a time series of 3 responses is used 
to extract 6 (i.e., 2 discrete and 4 continuous) parameters. In what follows we will use the 
notation 

Y, EE (yiiU), y2{U), ysiU))' 

and 

H^ = {^i{U,£,m,f3), ^i2{Uj,m,f3), ^■i{Uj,m., (3))' 

for the first three observed and theoretical moments respectively, at time point ti {i = 
1, • ■ • , n), where f3 = (vp, a, Ve, a)' . 

It is important to understand how the moment method is currently used. Theoretically, 
it can be shown (Aerts et al., 1992) that for a monopcriodic star, the time dependence of 
the moments takes the following form: 

/ii = fli sin(27ri'i + Ki), '1 
A*2 = bo + bi sin{2Tnyt + di) + b2 sm{ATrvt + 52), > (2) 

= ci sin(27ri^i + 7i) + C2 sin(47ri^t + 72) + C3 sin(67rzyi + 73), J 

where v is the oscillation frequency. The phases ai. Si, jj are constants, while the pos- 
itive amplitudes ai, hi, Cj depend on the parameters (£, to,/3). A discriminant r™(/3) is 
constructed to estimate these parameters by comparing the observed amplitudes with their 
theoretical counterparts: 

rr = I (/a, \ai - ail) + ^ (^/^^ vf^) + E </ \^'^ ~ ^'1 M ' (3) 

where the tilde denotes observed quantities, and where the weights / are introduced to 
incorporate the estimated standard errors Adi, Abi and Acj of the corresponding observed 
amplitudes: 

/ ^ VF-i / = i, h. = A, 

Aai' At, Aq' 
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(Aerts 1996). The form of F™ in © prevents the third moment 1/3 with its large values from 
dominating the first moment yi, but has the disadvantage that it cannot discriminate the 
sign of the mode number m. The parameters are estimated by searching for the minimum 
of in a rectangular grid in the parameter space. For the continuous parameters, it is 
hoped for that the grid is fine enough in order not to miss the global minimum. Finally, a 
table is produced with the top 5 or 6 best fitting {£,m,P) parameter sets. 

The best strategy to obtain the final estimate for (£, m) and f3 together with their 
uncertainties, is currently open to debate. Despite the usefulness of the table with the best 
parameter sets, there is no estimate of the uncertainties of the parameters obtained with 
the moment method, because of severe theoretical and computational complexity. This 
paper takes an important first step towards estimating the uncertainties of the continuous 
parameters f3. Estimating the uncertainties of the discrete parameters £ and m is an even 
more challenging problem, and will be left for future research. 

Even for a given {£, m) value, it is currently unknown how precise the continuous pa- 
rameters are estimated. For example, is the uncertainty in the inclination angle as small as 
5°, or is perhaps 30° a more typical value? Moreover, very often several {£,m) pairs give 
almost equally good fits. The question is raised as to how should we take this into account 
for our best estimate of /3 and its uncertainty? In what follows, we will try to answer these 
questions. 

4. New Statistical Approach 

We consider a new estimating method which produces both point and interval estimates. 

We first note that the three responses yi, j/2 and j/3 are dependent, and that their 
covariance matrix V is unknown. Formulating a statistical model for the noise on the 
moments is non-trivial as it would involve a model for both the instrumental and the 
atmospheric noise. In addition, we note that the relation between the coefficients oi, bi, 
Ci in (01, and the parameters (3 is non- linear thereby preventing the easy computation of 
a Jacobian matrix. Estimating /3 and its covariance matrix Cov[/3] using a simple variable 
transformation technique is therefore not possible. 

A first alternative is the least squares method. Although multi-response least-squares 
estimation has been used before to deal with correlated responses where the covariance 
matrix has to be estimated, Seber and Wild (1989) show that this technique should not 
be used if the covariance matrix V depends on the parameters /3, as is the case here. For 
example, the uncertainty in the first moment of a line profile (j/i) can be estimated with 
the second moment (2/2), and the latter depends on m, and /3. Or, with an astrophysical 
example, the faster the star rotates (larger Ve) the broader and flatter the line profile, and 
the less precision with which wc know the position or the first moment of the line profile. 
To avoid confusion, we stress that our argument in the example above is not that the 
uncertainty of the first moment yi depends on the uncertainty of the second moment y2 , 
but that the uncertainty of the first moment yi always depends on the second moment itself 
of which we know that it depends in turn on the parameters f3. 

We must therefore conclude that in the case of the moment method, minimizing the 
weighted sum of squares is not appropriate, regardless of how V is estimated, because it 
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will not yield a consistent estimate of the parameters f3 and it will reduce the efficiency of 
the estimator (Seber and Wild, 1989). 

The generalized estimating equations (GEE) methodology, as developed by Liang and 
Zeger (1986), is better suited for the purpose of the moment method. We recall that this 
method does not assume a particular joint probability density for the responses j/i, y2 
and 2/3, nor that they are i.i.d. The theory does not assume that the theoretical model is 
linear in its parameters, while the covariance matrix V of the responses does not need to 
be completely specified. The method does assume, however, that the different observations 
Yj (z = 1, • • • , n) are independent, that a working approximation of the covariance matrix of 
the responses is available, and that the expectation values -EiYi] = /i.,;(^, m, f3) {i = 1, ■ ■ ■ ,n) 
are correctly specified. 

We therefore use GEE to estimate the uncertainties of the continuous parameters f3. 
We recall that in the GEE method, the parameters are estimated by locating the root of 
the quasi-score function U(/3): 

N 

U(/3)^^D*.Wri.(Y,-M,), (4) 

where N is the size of the time series. The 3x4 matrix D = dfi/dfi^, and the 3x3 
symmetric matrix is a working approximation of the true covariance matrix of the 
quantities Y^: 

V, EE E[iY, - M,(/3))(Y. - tiMYl (5) 

where (3 are the true (but unknown) parameters. It can be shown (e.g., Liang and Zeger, 
1986; Zeger and Liang 1986; Diggle et al. 2002) that the root ^ is a consistent and asymp- 
totically normal estimate of the true /3, with sandwich covariance matrix 

Cov[^] =1-1 Ii Igi, (6) 

N 



where 



d(3 



and 

Ii Cov[U(/3)] = ^D* Wri V, Wr^ D, 



N 



The unknown covariance matrices in the expression for Ii are estimated by 

V, = (Y,-M,(;9)).(Y,-/x,(;9))*. 

The so-called sandwich estimator in Q is robust against misspecification of the covariance 
matrix of the responses. 

For the working approximation Wj for the covariance matrix , we suggest the follow- 
ing idea, where we estimate the uncertainty of the first three moments of the line profile 
with the higher moments, as is sometimes done with the moments of a probability distri- 
bution function. Consider the mirror image ({(j), v) ~ I —p{4>, v) of the spectral fine p{4>, v) 
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as a distribution function for the velocity v, and compute, for a given time U 



EC, 



E 



^ Cj [v] - Mr) ■ XI («| - Ms) 

] 3 

r • {iJ-r+s — l^r Ms)j 



where the sum over the index j runs over all the velocity points (pixels) of the spectral line, 
and where we define 




Here, we assume that the different observed points of the line profile are uncorrelated. The 
extra factor F appears because, contrary to probability distribution functions, line profiles 
are not normalized in area. Note that we use the higher theoretical moments and not the 
observational ones because, as mentioned before, the latter are often too noisy. It is difficult 
to assess the influence of the working approximation W on the final uncertainties on /3, but 
we refer to Diggle et al. (2002) where it is shown that the sandwich estimator © for the 
covariance matrix of /3 is quite robust against misspecification of V. 

Having derived an estimator for /3 and its uncertainty, given an (£, to) pair, we should 
take into account that we do not actually know the correct {£, m) values. If £ and to were 
continuous parameters, we would have a total of 6 continuous parameters for which we 
would have liked to compute a 6-dimensional confidence region. As I and to are discrete, 
however, it is notoriously hard to find an equivalent "confidence region" . We remind that 
it is current practice simply to take the (3 values of the best-fitting (^, to) pair with no error 
estimate at all. As a first alternative, we propose to "weight" each mode (£, to) with a 
lack-of-fit function. The best guess for both /3 and its uncertainty is then computed with 
a weighted mean over all relevant modes (£, m). The entire estimation procedure can be 
summarized as follows: 

(a) Specify a set of pairs of the degree H. and the azimuthal number to: {(^j, Wj)}. 

(b) For each of the pairs {£j,mj), solve the quasi-score equations and estimate the con- 
tinuous parameters f3j and their covariance matrix Cov[/3j]. 

(c) Compute for each of the modes {£j, rrij), the lack-of-fit parameter G| which indicates 

how well the theoretical moments /x(/3) fit the observed moments y: 



N 



Gi 



EE 



(d) The best estimate for the degree and the azimuthal number {£, to) is the {£j 



(7) 



that 



has the lowest lack-of-fit G^. The corresponding best estimate for the continuous 
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parameters (3 can be computed with 



and the corresponding covariance matrix is the sum of the intra-mode variance and 
the inter-mode variance: 

E Cov[/3;.] E ' ^j) ■ - GJ' 
Cov[P]='-^^ . (9) 

{iej,m.j)} {(^i.m,-)} 

For both practical and astrophysical reasons, only modes with a degree up to a certain 
limit (e.g. £ < A) are considered. 

In the following section this estimation procedure is applied to a dataset of the star 
HD181558. 



5. Application to HD181558 

HD181558 belongs to the class of the Slowly Pulsating B stars (SPBs). Although the star is 
multi-periodic (De Cat and Aerts 2002) it has a very dominant (in amplitude) first mode, 
which justifies a monoperiodic approximation. The amplitude of this mode is the largest 
ever observed for an SPB. The dataset used for this GEE application has already been 
shown in Figure 13 In what follows we always assume the theoretically predicted value 
K = 21. 

Our first goal was to estimate (3 for each mode {£, m) with < 4, by solving the non- 
linear quasi-score equations. It turned out, however, that this was not just a technical detail 
of the procedure, but was in fact a major issue. 

First, it turned out that the quasi-score function U(/3) is computationally slow to eval- 
uate, with one evaluation requiring 18 time series evaluations: 6 for the moments fJ-i-fJ-e for 
the working approximation W, and 12 for the moments /X1-/X3 for different parameters (3 to 
numerically compute (with forward differences) the derivatives in D. For this reason, prior 
to using we first determined a good initial guess for ^ for the local search routine, using 
a rough scan of the 4D parameters space for each mode {i, m) with a computationally less 
expensive lack-of-fit function g{(3): 



\ 



1 ^ 

j^Y. \yd{u)-iid{u,/3)\. (10) 



The construction with the d*'' root and the division by d simply prevents the higher order 
moments from numerically dominating the lower order moments. The sampling of the 
parameter space was done probabilistically and non- uniformly. For each parameter (3i, a 
physical range was determined and this range was subdivided into intervals. After each 
set of 10000 sampled points, each interval of each parameter (Ji was assigned a sampling 
probability according to the lowest g{0) value recorded up to then, with the f3i component 



10 J. De Ridden G. Molenberghs and C. Aerts 

Table 2. For each {£, m) pair, 
the 4D parameter space was 
scanned with the lack-of-fit func- 
tion g defined by equation 
and with the dataset of the star 
HD181558 shown in Figure 
The minimum gmin of the g{/3) 
values for each of the {£, m) pairs 
is given. 



m 










1 


2 


3 


4 


+4 








11.9 


+3 






7.52 


11.0 


+2 




6.57 


6.71 


11.3 


+1 


4.72 


4.74 


5.85 


10.8 





6.37 


6.37 


6.79 


11.6 


-1 


4.79 


6.57 


7.05 


10.7 


-2 




4.68 


6.86 


10.5 


-3 






6.92 


11.3 


-4 








11.3 



in the corresponding interval. The sum of probabihties over aU intervals of a parameter f3i 
was set to one. For each {l,m) pair a total of 200,000 points was sampled. This procedure 
was set up to sample the more promising regions of the parameter space. 

In Table[5]we give for the star HD181558 the lowest g{l3) value recorded for each mode 
{£,m). The {i,ni) — (0,0) pair can be excluded on astrophysical grounds because such 
modes do not occur in SPBs. As can be seen, no mode {£,m) stands out, but there are 
several candidate modes that describe the data well. Our final estimate of /3 should take 
into account this mode uncertainty. We also note that Table|21is not symmetric with respect 
to the sign of m. The moments indeed behave differently when the pulsational wave goes in 
the same direction as the rotation then when the wave goes in the opposite direction. This 
sensitivity to the sign of m was lost in the old approach (outlined in Section |3J) where one 
only uses the absolute value of the amplitudes. 

The 24 scans of a 4D parameter space with 200,000 points each, was a rather time 
consuming but necessary task to find suitable initial guesses for (3 for the local search 
algorithm. We implemented two derivative-free methods: the conjugate-direction (Powell's) 
method (see, e.g., Press et al., 1992, p. 420) and the Torczon (1989) simplex method. The 
former of which, having the best performance, was used to locate the root /3 of U for all 
modes {£,m). We found that these methods had much stabler performance than quasi- 
Newton methods, such as Newton-Raphsons, Fisher scoring, or variations to this theme. 

Even with the conjugate-direction method, the algorithm did not always converge. The 
reason, as it turns out, is that the quasi-score functions have "false" zeros, for example there 
are cases where the components of U approach zero for cr — > oo. Quite often, the algorithm 
converged to a point outside the physically relevant range of the parameters, even when 
several different initial guesses for ^ were tried. Although they did not occur for our dataset 
of the star HD181558, we should mention two other possible causes of numerical difficulties. 
First, it may be possible that the working approximation W is not invertible, for example 
if f3 approaches zero. Second, the matrix 1q may not be invertible, and hence no covariance 
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matrix can be computed. This occurs, for example, for a — > 0° because Ve appears only 
in Ve sin a in the equations, so that the third row and the third column of lo are zero. We 
stress, however, that the latter example is a problem of intrinsic non-idcntifiability and is 
not specific for the GEE approach. One simply cannot derive the rotational velocity if the 
star is looked pole-on. 

Making detailed ID slices of the 4D function ||U|| too time consuming, but we record the 
minimal lack-of-fit values gmm in each of the intervals of each parameter (disregarding the 
values of the other parameters). Figure O shows typical examples. Although the function 




Fig. 3. Representative examples of the minimal lack-of-fit value gmin for each sample interval of a 
parameter, for the star HD1 81 558. We remark that although the size of the intervals for the parameter 
Vp is fixed, the relevant range of Vp depends on the mode numbers {£, m). 

g{(3) need not have exactly the same behaviour as the function ||U(/3)|| (the difference 
is similar to the well-known difference between Li-norm and L2-iiorm minimization), we 
assume that the functions share many features. We observe that the minimum in the upper 
left panel of Figure for the well-fitting mode {£,m) = (2,-2) is quite localized. This is 
much in contrast with the almost flat surface in the lower left panel for the badly-fitting 
mode {£, m) = (4, —3). Intuitively, one can expect that the equivalent for the case of the ||U|| 
function hampers the iterations towards the minimum, and that this increases the chance 
of wandering from the physically relevant part of the parameter space. This is exactly what 
happened for this mode. More generally, we observe strong correlation between how well 
a mode fits the data (with gmin as lack-of-fit value) and the chances that the root finding 
algorithm does not converge. The lower right panel shows two minima in the plot of the 
inclination angle a, since we consider a full 360° range. While over such a range symmetry 
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Table 3. Roots of the quasi-score functions for those modes where there was convergence 
in the physically relevant part of the parameter space, for the star HD1 81 558. The values 
between brackets are the standard errors obtained with the sandwich estimator is 
the lack-of-fit value of the mode as defined by Eq. Q. Vp, a and Ve are expressed in km/s, 
and the inclination angle a in degrees. 



il,m) 


lluiLi„ 




Vp 


a 




a 


(1,0) 


0.15 


2.7 


2 (1) 


9.0 (0.9) 


13 (19) 


320 (48) 


(1,1) 


5.1 10-2* 


0.63 


2.01 (0.08) 


6.3 (0.2) 


15.2 (0.7) 


117 (2) 


(1,-1) 


1.1 10-5 


1.1 


4.0 (0.2) 


4.2 (0.6) 


25 (2) 


336 (1) 


(2,0) 


0.10 


3.1 


0.9 (0.4) 


7(2) 


30 (24) 


331 (13) 


(2,1) 


1.7 


0.61 


1.7 (0.1) 


3.8 (0.8) 


16 (1) 


71 (1) 


(2,2) 


0.014 


2.4 


1 (1) 


10 (2) 


0.4 (29) 


270 (360) 


(2,-2) 


0.0032 


0.72 


1.62 (0.06) 


4.3 (0.4) 


17.6 (0.7) 


129 (1) 


(3,1) 


1.1 10-25 


2.3 


1.00 (0.04) 


6.8 (0.7) 


18 (2) 


145 (7) 


(3,2) 


0.40 


3.1 


1.1 (0.2) 


6(3) 


17 (10) 


49 (23) 


(3,-1) 


3.5 


7.0 


1.3 (0.3) 


3 (12) 


7(12) 


189 (5) 


(4,0) 


2.7 


11 


0.4 (0.3) 


8 (12) 


49 (188) 


19 (96) 


(4,-4) 


0.047 


8.7 


0.6 (3) 


7 (24) 


20 (87) 


295 (30) 



relations exist, these depend on i and m and taking them into account to limit the range of 
a would entail a lot of bookkeeping that can elegantly be avoided by simply considering the 
entire range. The feature that we quite often do not seem to find the root of U does not 
necessarily contradict the theory outlined in Section^ We solve for the root of the observed 
U function because we know that E[\J] = 0. However, the latter is only true if the model 
is correctly specified, i.e. if E{Y] = fi{£,ni,(3). Therefore, theoretically, the existence of a 
root in the 4D parameter space of the continuous parameters (3 cannot be guaranteed for a 
"wrong" {£, m) pair, and this is exactly what we observe for badly fitting modes. For this 
reason we interpreted a non-convergence (after repeatedly trying) as an indication that the 
candidate mode should be disregarded. 

In Table 01 we list the roots of the quasi-score function for those modes for which there 
was convergence in the physically relevant part of the parameter space. The closeness of ||U|| 
to zero varies from mode to mode. In cases where this value fails to be small, we checked 
this is not due to premature convergence, since restarting the algorithm at the point where 
it stopped, did not further decrease ||U||. One possible explanation might be that, for some 
of the solutions, the algorithm has converged to a local minimum, such as for the modes 
(£,m) = (3,-1) and (4,0). 

In Figure ^ we show fits for the three best fitting modes, with the function (see 
Eq. d) as a lack-of-fit. As mentioned before, there is not just one, but several modes 
that can fit the observed data quite well. The lack-of-fit values in Table El provide us 
with an indication of the relative merits of wave number choice {£,m). Of course, at 
this point we lack knowledge about the reference distribution of these values, unlike in 
classical fit statistics (e.g., likelihood-ratio based). However, similar instances exist in both a 
frequentist (e.g., Akaike Information Criterion) and a Bayesian context (e.g., Bayes factors). 
Nevertheless, we assert that these numbers, especially when supported by careful graphical 
inspection, are useful to narrow down substantially our uncertainty about the wavenumbers, 
in spite of an intrinsically complicated modelling endeavor. To this end, the last column 
in Figure 01 displays mode {£,ni) = (3, 1) with a substantially worse fit than the one in the 
first three columns of the same figure. 
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(l,m) = (1.1) (l.m) - (2.1) (l.m) - (2.-2) (l,m) - (3.1) 



4000 - . 




0.5 1 0.5 1 0.5 1 0.5 1 

phase phase phase phase 



Fig. 4. Theoretical models (solid lines) of the observed moments (bullets) for the three best fit- 
ting modes {£,m) = (1,1), {£,m) = (2,1), and {e,m) = (2,-2), plus the poorer fitting mode 
{e,m) = (3,1), with as a lack-of-fit function. The theoretical models were obtained with the 
model parameters obtained with the GEE method. The first, the second and the third row are for the 
first moment yi (km/s), the second moment y2 (km^/s^) and the third moment ys (Wm^/s^) respec- 
tively. The moments are shown as a function of the phase. Note that the models of the different 
promising modes differ mainly for the second moment. 



We used the modes in Table 13 to compute the weighted mean /3 and its standard error, 
with Eqs. ||HJ and 0. The results, including the intra- mode and inter-mode variance, are 
given in Table 01 In the specific case of HD181558, one could argue that the modes with 
£ = 3 and £ = 4 can be disregarded on astrophysical grounds. The reason is that these 
modes would require a very large oscillation amplitude at the surface of the star to cause 
the large observed amplitude of the first moment yi. For this reason, we also computed (3 
with the £ = 1 and £ = 2 modes of Table |31 only. The results are listed in Table [S] With 
Tables 01 and [SI we achieve the goal of this application of the revised version of the moment 
method: we have obtained a best guess for the continuous parameters and their standard 
errors, where we took into account the mode uncertainty. An important result is that the 
uncertainties of the parameters can be large, in fact larger than we expected. Especially 
the rotational velocity cannot be estimated precisely. The large inter-mode uncertainty 
of the inclination angle a is not surprising since the inclination angle is known to be largely 
dependent on the mode numbers {£,m). We note that the values for the weighted means 
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Table 4. The weighted mean over all 1 2 modes in Table|31 computed with 
Eqs. js) and |9|. The values mentioned between brackets are standard 
errors. The intra-mode variance and inter-mode variance are computed 
with respectively the first and the second term of Vp, a, Ve are ex- 
pressed in km/s, and the inclination angle a is expressed in degrees. 



ft 


Weighted Mean 


Intra-mode Variance 


Inter-mode Variance 


Vp 


1.8 (1.0) 


0.33 


0.74 


a 


5.5 (4.1) 


14 


3.1 


Ve 


17 (26) 


612 


46 


a 


164 (132) 


7300 


10079 



Table 5. The same information as in Table|4]is shown, except that the 
mean is computed over those 7 modes in Table|31with £ = 1 and £ ^2. 



ft 


Weighted Mean 


Intra-mode Variance 


Inter-mode Variance 


Vp 


2.0 (1.0) 


0.19 


0.71 


a 


5.3 (2.0) 


0.82 


3.1 


Ve 


16 (12) 


100 


37 


a 


170 (137) 


8342 


10450 



do not change much by excluding the £ — 3 and £ — A modes. The reason is that the latter 
modes have a lower weight anyway, as can be seen from the values in Table |3I 

The fact that the standard errors for the continuous parameters are large turns out not to 
be specific for the star HD181558. We applied our method to several artificial datasets, and 
we obtained similar results. Hence, we conclude that estimates of the continuous parameters 
generally can be quite uncertain. 

6. Summary and Conclusions 

We made a first but arguably important step to develop a statistical formalism for the 
moment method. In this first stage we aimed to incorporate estimates of the uncertainties 
of the continuous parameters. Because of the many difficulties to overcome, this was never 
done for the moment method, nor for any other mode identification technique. 

We found that, in the specific case of the moment method, the method of least-squares 
does not give consistent estimates of the continuous parameters and we resort to the GEE 
method (Liang and Zeger 1986). This method requires a working approximation of the 
covariance matrix of the 3 responses, based on the higher theoretical moments. Note that 
the higher moments, known to be imprecise, are not used in the actual model. An important 
source of uncertainty is the fact that often not just one but several candidate modes can 
describe the data. We set up a separate procedure to weight each mode and to compute a 
weighted mean over all modes of the parameter vector and its uncertainty. To compute the 
latter we introduced the intra-mode and the inter-mode uncertainty. 

Subsequently, we apphed our procedure to the SPB star HD181558, from which we 
learned the strong and the weak points of our method. We found out that solving the 
estimating equations is a computationally demanding and tedious task: convergence of 
the algorithm was not evident, despite the fact that we experimented with several robust 
local root finding methods, of which we selected the method of conjugate directions as 
the most efficient one. On the other hand, we also proposed a new lack-of-fit function to 
scan the parameter space to obtain good initial guesses for the local search method. This 
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lack-of-fit function proved to be very useful on its own, as it works at least as well as the 
old discriminant (O and allows in addition to discriminate between positive and negative 
azimuthal numbers which was one of the shortcomings of the previous discriminant. Our 
strategy to scan the parameter space also proved that there are several modes that can 
explain the dataset of HD181558. Only taking into account the very best fitting one, would 
therefore not be useful. 

This is why we retained 12 modes as candidate modes for which an estimate of the 
continuous parameters f3 can be computed, and used these estimates to obtain a best guess 
(3 for the continuous parameters plus their uncertainties, taking into account the mode 
uncertainty. Doing so, we discovered that the parameter uncertainties can be large, a result 
which was moreover confirmed in the case of artificial datasets. 

Prior to our study, such large uncertainties were not anticipated. On the contrary, 
one sometimes assumed the uncertainties to be quite small in order to be able to apply a 
two-stage approach in the multiperiodic case. In such an approach the inclination angle a 
and the rotational velocity Ve are determined with the dominant mode, and subsequently 
fixed while determining the mode parameters of the other modes, to have the dimension 
of the parameter space reduced. Our new results show that such an approach can be 
very dangerous: in the case of HD181558 it can hardly be justified because of the large 
uncertainty on a. 

The method we outlined in this paper is the very first attempt to develop a statistical 
formalism for the moment method. Even though there is undoubtedly room for additional 
work before our proposed method can be deemed widely applicable, we conclude it makes 
an important first step in our understanding of the uncertainties and usefulness of the 
continuous parameters, estimated with the moment method. It furthermore underscores 
that some conclusions reached in the past need to be revisited. 

We conclude with several possible future improvements. First, it may be worth to in- 
vestigate alternative parameterizations /3' which, while mathematically equivalent to the 
original one, may improve upon the convergence and robustness properties of the algo- 
rithm. We already experimented, for example, with using VeSina instead of Ve, and with 
using Vp sin a and Vp cos a instead of Vp and a. Second, it would be useful to extend the 
formalism to include the uncertainty on K. Third, it might also be interesting to use several 
spectral lines at the same time to improve the statistics. Including multiple modes might 
also improve the convergence properties. Although this would imply 3 more parameters 
(Z2, r7i2, Wp,2) per mode, multiple modes would set more stringent restrictions on the incli- 
nation angle a which has a significant impact on the estimation of the other parameters. 
To further validate all of the above, more simulations are necessary. In addition, such 
simulations can also clarify what impact the number and the signal to noise ratio of the 
observational spectral lines has on the performance of the algorithm. 
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